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Abstract 

The interaction of submerged manuevering vehicles with the free surface is a topic of im- 
portance in ocean engineering. Present methods for estimating the forces and motions due 
to wave-body interactions are limited in their ability to efficiently compute forces and mo- 
tions of interest when free surface slopes are steep or body motions are large. Methods are 
required which contribute to overcoming these difficulties. 

This thesis presents a method for analyzing the nonlinear interaction of non-breaking 
waves with submerged bodies of arbitrary geometries undergoing arbitrary motions. The 
method couples a high order spectral representation of the free surface profile and free 
surface potential with a boundary element representation of the body to effect a time domain 
solution for specified initial-boundary value problems. The spectral representation of free 
surface quantities facilitates the use of fast transform techniques for rapidly computing 
free surface quantities. The boundary element repesentation of the body enables bodies of 
arbitrary geometry to be modelled. Fully nonlinear free surface boundary conditions axe 
used as time evolution equations for the free surface profile and potential. 

The effectiveness of the method is demonstrated through the solution of three classes 
of two-dimensional problems: (1) diffraction of incident waves by a stationary cylinder, (2) 
wave radiation by a cylinder undergoing forced oscillations, and (3) combined diffraction 
and radiation by a submerged neutrally buoyant cylinder free to respond to incident waves. 
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Chapter 1 



Introduction and Background 



This thesis presents a method for solving wave body interaction problems. The method 
developed is appropriate to the analysis of the nonlinear interaction of non-breaking waves 
with fully submerged bodies of arbitrary geometry undergoing arbitrary motions. The 
method is suitable for the study of three or two dimensional problems. The efficacy of the 
method is demonstrated through the study of three classes of two dimensional problems: 



• The diffraction of an incident Stokes’ wave by a fully submerged cylinder of arbitrary 
shape; 

• The wave radiation by a fully submerged circular cylinder undergoing large amplitude 
forced oscillatory motion; and, 

• The response of a neutrally buoyant circular cylinder to an incident Stokes’ wave 
(combined radiation and diffraction). 

The purpose of this chapter is to explain the motivation, background, and format for 
this thesis. 
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1.1 Motivation 



This thesis is motivated by the desire to efficiently estimate and understand the interaction 
of fully submerged bodies with the free surface. A time domain solution for the forces 
resulting from the interaction of the body with waves is necessary to perform maneuver- 
ing simulations of bodies operating near the free surface. Present methods for estimating 
forces due to wave-body interactions limit the computing speed of maneuvering simulation 
calculations or include linearizing assumptions which limit the ability to predict forces and 
motions of interest when free surface slopes are steep or unsteady body motions are large. 
This thesis presents a method which overcomes some of these shortcomings. 



1.2 Background 

The Physical Problem 

The complete physical problem for a body operating in proximity to a free surface is 
dominated by inertial, viscous, and gravitational forces [33]. Other mechanisms, such as 
surface tension, exist but are insignificant. A mathematical representation of the complete 
physical problem amenable to efficient solution does not exist. The problem includes the 
possibility of separated flows and breaking waves, phenomena difficult to analyze. Numerical 
methods enable a range of problems that include some of these phenomema to be solved, 
however the computing time required by some methods to solve these problems for arbitrary 
bodies is prohibitively excessive. For example, the solution of the Navier-Stokes equations 
for the flow around a submarine in an unbounded fluid for six body lengths of travel requires 
as much as 100 hours of CPU time on an IBM Model 590 workstation with 512 MB of RAM 
[27]. Numerical methods are needed which can characterize the relevant fluid and body 
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behavior with improved computational efficiency. 

The relative importance of the principal force mechanisms varies as the parameters 
characterizing the problem change. Dimensional analysis is used to define the relationships 
between the force mechanisms and identify ranges of problem parameters over which one or 
more of the force mechanisms may be ignored. For a typical marine vessel moving beneath 
the free surface the ratio of inertial to viscous forces is large outside the boundary layer. 
The fluid can be treated as inviscid and the governing mathematical representation of the 
problem appropriately simplified [33]. If the problem is further restricted to considering 
only non-breaking waves, a wide range of practical problems can be solved. In the following 
section, the existing methods for solving wave-body interaction problems for inviscid fluids 
under non-breaking waves are briefly discussed. 

Solution Methods 

The simplest solution method is that proposed by Faltinsen [9]. This method approxi- 
mates the hydrodynamic force as the sum of (1) the Froude-Krylov Force (Ff*;), the force 
due to the unsteady pressure field associated with the undisturbed incident wave potential; 
(2) a force proportional to the product of the body’s added mass and the acceleration vector 
of the undisturbed wave field; and (3) the added mass force, the force resulting from the 
body’s acceleration. This approach has been used [5] for simulating the forces associated 
with undersea vehicles operating near the free surface. The shortcomings of this type of 
approach are discussed by Newman [32] and Milgram [29]. This approach can only be 
accurate if 

• The Froude Numbers based on submergence depth, F//, and body length, Fi are 
small: 
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F h = -¥= <1 and F l = -^= < 1 

VgH vgL 



( 1 . 1 ) 



where is the body velocity, g is the gravitational constant, H is the body’s mean 
submergence, and L is the body’s length; and, 

• The body is small compared to the wavelength. 

Because of these restrictions, the method is not suitable for many combinations of marine 
vessels and operating scenarios of interest. 

Numerical solution methods for wave-body interaction problems exist in both the fre- 
quency and time domains. The application of the principle of linear superposition to ship 
motions in irregular seas was motivated by the work of John [15] and St. Denis and Pierson 
[44]. The former, as noted by Korsmeyer [19], enabled the problem to be decomposed into 
separate radiation and diffraction problems, and the latter, as noted by Salvesen, et al. 
[40], enabled ship responses in irrregular seas to be evaluated as the sum of the responses at 
all frequencies of regular waves. Jointly these works motivated the development of a wide 
body of linearized ship motion analysis methods in the frequency domain using strip theory 
and slender body theory to characterize the wave-body interaction. The work of Salvesen, 
Tuck, and Faltinsen [40] is a thorough application of strip theory to ship motion analysis. 
Ogilvie [36] provides a thorough development of slender body theory and its application to 
ship motions. 

Improvements in digital processing capabilities enabled the development of three dimen- 
sional analysis methods for arbitrary wave induced ship motions [2], Solution methods have 
evolved systematically, from the solution to the linearized problem for stationary bodies or 
bodies undergoing forced motion to the solution for bodies undergoing nonlinear unsteady 



16 



motions with forward speed. Both frequency and time domain solution methods have de- 
veloped for such problems. The solution to the specified boundary value or initial value 
problems is found by either volume methods (finite difference and finite element methods) 
or boundary integral equation methods (boundary element or “panel” methods) [7]. The 
present preferred methods in ship hydrodynamics are the boundary integral equation meth- 
ods (BIEM). Kring [20] summarizes the development of three dimensional panel methods, 
with emphasis in his summary on the distinction between methods which use the transient 
free surface Green function ( Korsmeyer [19] and Bingham [2]) and those which use the 
Rankine source (Nakos [31] and Kring [20]). Time domain simulations of wave-induced mo- 
tions are of immediate interest for their applicability to vehicle motion control and control 
system design. 

Dommermuth [7] and Liu [22] developed a high-order spectral method (HOSM) for the 
efficient time-domain analysis of non-linear wave-body interactions. The high-order spectral 
method, through the use of Fast Fourier Transform (FFT) techniques, greatly improves the 
computational efficiency of the numerical analysis of wave-body interactions. 

The solution method developed in this thesis couples the boundary integral representa- 
tion of the body with a spectral representation of free surface quantities. To demonstrate 
the efficacy of the method, the three classes of problem identifed previously, radiation, 
diffraction, and combined radiation and diffraction, are studied for two dimensional bodies. 
The solutions are compared with solutions generated by different methods for similar two 
dimensional problems. These other methods will be briefly summarized here. A similar 
summary of some of these methods and related methods can be found in Dommermuth [7] 
and Liu [22]. More detailed discussions of the methods and the solutions generated will be 
discussed in subsequent chapters as appropriate. 
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Frequency domain solutions to wave diffraction and radiation by circular cylinders were 
developed by Ursell [47] and Ogilvie [35]. Ursell used a multipole expansion of the pulsating 
source potential and a series expansion, in powers of the wave slope, for the coefficients 
in the multipole expansion, to derive a system of equations to uniquely determine the 
linear solution to the interaction between waves and a fully submerged stationary circular 
cylinder. Ogilvie extended Ursell’s work to compute the first-order unsteady force and the 
second-order steady force for the three problems of interest in this thesis. Wu [53] extended 
Ogilvie’s work using the multipole expansion and the free surface boundary condition for the 
second order potential to demonstrate that to second order, a circular cylinder undergoing 
a circular orbit generates waves in one direction only. Ogilvie previously demonstrated this 
result at first order. Using a linear free surface boundary condition and the exact body 
boundary condition, Wu [52] found that a circular cylinder undergoing large- amplitude 
circular motion can generate waves travelling in both directions though the waves travelling 
in the “upstream” direction were at most a third order quantity. 

Chapman [4] analyzed large amplitude transient motions of two-dimensional bodies with 
a linearized free surface boundary condition and exact body boundary conditions. Chapman 
represented the wave field as a finite sum of eigenfunctions and the body as a distribution of 
sources and negative image sources reflected about the mean free surface. For Chapman’s 
method the periodic boundary conditions required by the spectral representation of free 
surface quantities are ignored in the representation of the body potential, the harmonic 
components used to represent the free surface quantities are nonuniform and can not be 
determined by FFT, and the method is applicable to flows described by linear free surface 
boundary conditions. 
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Grue and Palm [11] studied the reflection of waves by submerged cylinders of circular and 
elliptic cross section using boundary integral equations. In their formulation the solution 
is expressed as a distribution of vortices on the body. Grue and Palm used the linear free 
surface boundary conditions in their formulation and confirmed the results of Ursell and 
Ogilvie (no reflection from a circular cylinder) and computed the reflection power for the 
elliptic contour. Their analyses of cylinders of elliptic cross section are restricted to deeply 
submerged cylinders. 

Stansby and Slaouti [45] developed a time domain solution method for evaluating non- 
linear wave-body interactions with the free surface represented as a periodic vortex sheet 
and the body as a distribution of periodic sources. The focus of their work was on the first 
order oscillatory forces and wave motions only and they demonstrated reasonable agreement 
with the first order solution of Ogilvie, in both the amplitude and phase of the oscillatory 
forces. 

A numerical solution to the second-order wave diffraction problem in the frequency 
domain was provided by Vada [48]. Vada derived the boundary conditions for the first and 
second order diffraction potentials using perturbation expansions of the boundary conditions 
and solved the first and second order problems using a boundary element method. Vada used 
the frequency domain free surface Green function provided by Wehausen and Laitone [50] 
for the first order problem and a suitable modification of it for the second order problem. 
The results generated by this method compare favorably with the results of Ogilvie and 
experimental results [3]. Vada did not solve for the time-independent part of the second 
order potential and restricted his study to stationary bodies. 

The Mixed Eulerian-Lagrangian (MEL) method has been applied to the study of non- 
linear wave body interactions by a number of authors. Vinje and Brevig [49] and Cointe [6] 
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studied the interaction of cylinders with non-linear free surface boundary conditions. Vinje 
and Brevig demonstrated the suitability of the approach to the study of wave-body inter- 
actions. Cointe studied the diffraction and radiation problems in detail and demonstrated 
the ability of non-linear solution methods to capture behavior not predicted by linear and 
second order theories. However, the Lagrangian representation of free surface quantities is 
not amenable to solution by fast transform techniques. 

Vinje and Brevig’s work, which used Cauchy’s integral theorem to represent the velocity 
potential and stream function , was used by Jagannathan [14] to study radiation conditions 
for the MEL method. In the process, Jagannathan produced time simulations for the forces 
on a submerged circular cylinder. The force results vary significantly from theortical and 
experimental results for wave diffraction but show better agreement for forced motions. Ad- 
ditionally, the difference between Jagannathan’s results and experimental results increases 
as submergence decreases. 

Silva and Peregrine [43] also used the MEL approach to study nonlinear wave-body 
interaction problems for submerged fixed bodies and submerged bodies undergoing forced 
oscillatory motions. Silva and Peregrine demonstrated the influence of finite depth on the 
forces and free surface profiles associated with forced body motion. Silva and Peregrine’s 
approach is not extendable to three dimensions. 

The high-order spectral methods (HOSM) of Dommermuth and Liu have been used to 
study nonlinear wave-body interactions in both three and two dimensions. By performing 
computations at orders higher than second order, Dommermuth and Liu were able to com- 
pute a steady horizontal drift force acting on a fixed submerged circular cylinder. This force 
is not predicted by second order theories but has been measured experimentally [30] and 
computed by fully nonlinear methods [6]. The HOSM uses the FFT to rapidly compute the 
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body and free surface potential. The HOSM requires that the body be capable of being 
efficiently represented by series of orthogonal functions. Additionally, the HOSM represents 
the body potential as a Taylor series expansion about the mean body position. Therefore, 
the method can not be used to model a body undergoing large amplitude unsteady motions. 

1.3 Thesis Format 

This thesis consists of seven chapters. The mathematical formulation of the combined 
high-order spectral boundary integral equation method is provided as Chapter 2. The 
numerical implementation of the method is provided in Chapter 3. The results for the radi- 
ation, diffraction, and combined problems are provided in Chapters 4, 5, and 6 respectively. 
Chapter 7 provides conclusions and recommendations for future work. 



21 



Chapter 2 



Approach 



This chapter presents the method developed for evaluating the interaction of a fully sub- 
merged body with the free surface. The problem statement, the solution formulation, and 
issues associated with the formulation are presented. 

2.1 Introduction 

The motivation for the mathematical formulation to follow is the solution of the equation 
of motion for a submerged body subject to both deterministic and stochastic forcing. The 
deterministic forcing is the result of control surface motions and operation of propulsive 
devices. The stochastic forcing is due to the interaction of the body and its appendages 
with waves. The problem formulation presented here provides an efficient method for com- 
puting the forces and moments on a fully submerged body of arbitrary geometry due to its 
interaction with waves. 

The formulation provides a mathematical model for the physical problem represented 
by Figure 2-1. 
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Figure 2-1: Coordinate system definition 



2.2 The Initial-Boundary Value Problem 

Field Equation 

For a large portion of the range of Froude numbers and Reynold’s numbers at which marine 
vessels operate while in proximity to the free surface, viscous effects are confined to a 
thin layer of fluid near the vessels’ surfaces. Prandtl’s boundary layer theory motivates 
and provides justification for neglecting viscous effects in the study of selected problems 
in hydrodynamics [33]. Chaplin [3] demonstrated experimentally that in the absence of 
viscous fluid effects significant non-linear wave forces exist. 

Wave-body interactions are dominated by fluid flows which can idealized as inviscid and 
irrotational. Therefore a portion of the non-linear aspects of wave-body interactions can be 
studied using potential flow theory. 
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For an incompressible fluid of constant density the equation for conservation of mass is 



the continuity equation 



V-V = 0 (2.1) 

For an ideal fluid initially at rest and subject to conservative forces the fluid rotation rate 
is zero and 



V x V = 0 



( 2 . 2 ) 



The velocity vector for an ideal irrotational flow can be defined as the gradient of a 
scalar potential 



V = 



(2.3) 



Equations (2.1) and (2.3) are combined to yield Laplace’s equation, the governing partial 
differential equation for the velocity potential: 



V 2 <F = 0 



(2.4) 
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Boundary Conditions 



The following boundary conditions are imposed on the physical problem represented by 
Figure 2-1: 



on 


on S Q 


(2.5) 


£>$ r . . 
— = U -n 
on 


on S B (t) 


(2.6) 


^ lx=0 ^ 1 x=2L 




(2.7) 


v* u„= v* |„ !t 




(2.8) 


O 

II 

1 

Cl IS 


on y = r] 


(2.9) 


d(b 1 

Pa) = + 2 V< ^ • V<£ + gy = 0 


on y = 7] 


(2.10) 



where i](x, z, t) is the shape of the free surface, U is the body velocity, and 2 L is the domain 
length. 

In the absence of a bottom: 

V$ -> 0 as y —> — oo (2-11) 

The specified form of the kinematic and dynamic free surface boundary conditions (equa- 
tions (2.9) and (2.10) respectively) determines the type of free surface flows represented. 
An Eulerian form of the free surface boundary conditions capable of describing non-linear 
single-valued free surface flows is used in this thesis. The chosen form was developed by 
Zakharov [55] for the purpose of evaluating the stability of dispersive surface waves. The 
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“Zakharov Equations” express the free surface boundary conditions as functions of the 
canonical variables <$ 5 (x,z,£) = $(£, 77 , z, t) and 77 , and the vertical velocity on the free 
surface, 3> y (x, 77, z, t): 



Vt = (1 +vl)$y(x,v,t) + $x T lx (2.12) 

<*>f = -m - \(* s x ? + 1(1 + T)l)e y {x, 7 7, t) (2.13) 

Spatial periodicity (equations (2.7) and (2.8)) is specified. The significance of this 
condition is explained in chapter 3. 



Initial Conditions 

The initial conditions for the the shape of the free surface, 77 ( 2 ;, z,0), and the free surface 
potential, <J> 5 (x,z,0), are specified as those associated with Stokes’ waves. 

2.3 Solution Method 

The Zakharov forms of the kinematic and dynamic free surface boundary conditions are 
integrated in time to update the free surface shape and free surface potential using the 
specified initial conditions and the solution of the boundary value problem. A combined 
high-order spectral method (HOSM) and boundary integral equation method (BIEM) is 
used to solve the specified initial-boundary-value problem. The Zakharov equations were 
used in the HOSM developed by Dommermuth and Yue [ 8 ] and West, et. al. [51] to 
solve non-linear wave- wave interaction problems and by Liu, Dommermuth, and Yue [23] to 
evaluate non-linear wave-body interactions. Olmez [37] used the Zakharov equations and 
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a boundary integral equation method to analyze non-linear wave- wave interaction problems. 



The high order spectral method (HOSM) and (BIEM) are summarized in the following 
sections. 



High-Order Spectral Method 

The HOSM summarized here is that of Liu, Dommermuth, and Yue [23] for solving nonlinear 
wave-body interaction problems. In the HOSM, and in the combined HOSM-BIEM solution 
method developed herein, the velocity potential is expanded in a perturbation series up to 
a specified order M : 



M 

$ {m) (x,y,z,t) = Y (2.14) 

771 = 1 

( denotes a quantity of order (e)( m ) where e is a measure of wave steepness, kA. 
Each perturbation potential, y, z, £), is expanded in a Taylor series about y = 0 

such that: 



$ S {x,Z,t) = $(x,7),z,t) 



M M—m 



E E 






(2.15) 



For the specified initial conditions, equation (2.15) establishes a sequence of Dirichlet 
boundary conditions for each perturbation potential, 0, z, t): 



$^(x, 0,z, t) = $ s (x,z,t) 



(2.16) 
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m = 2, 3, . . . , M 



(2.17) 



m-1 I pA 

$W(x,0,z,t) =- p l^^-0( Xf0 , 2 ,t) , 

The HOSM, as developed for wave-wave interactions, represents the velocity potential 
at each order as a Fourier series. The amplitude of each Fourier mode is determined, at 
each order m, by the Dirichlet conditions expressed above. In the HOSM developed for 
solving wave-body interaction problems, the velocity potential at each order is expressed as 
the sum of expansions of free surface and body basis functions [23]: 

oo oo 

(2.18) 

71=0 71=0 

where are the Fourier modal amplitudes of dipoles distributed on the free surface and 
a n (t) are the Fourier modal amplitudes of sources distributed on the body. 

The basis functions, typn and satisfy the Laplace equation, spatial periodicity 

requirements, and the bottom boundary conditions. For two dimensional problems, for 
example, the basis functions are given as Fourier integrals involving the two-dimensional 
periodic source potential. A derivation and explanation of the two-dimensional periodic 
source potential are available in Newman [34]. 

Boundary Integral Equation Method 

A boundary value problem governed by the Laplace equation can be recast as a boundary 
integral equation for the unknown velocity potential [13]. A review of the boundary integral 
equation method is available in Hunt [13]. 

The boundary integral equation method (BIEM) uses Green’s theorem to derive an 
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integral equation for the velocity potential in terms of distributions of singularities on the 
boundary of a simply connected fluid domain. The velocity potential is expressed as a 
function of singularity distributions on the fluid boundary by applying the second form of 
Green’s Theorem [12] 




dS = 0 



(2.19) 



where the Green function G and the velocity potential $ are solutions of the Laplace 
equation in the fluid domain. 

If the Green function is the potential of a source ( G = logr), equation (2.19) can be 
replaced by [13], 

/ L * (?)2 Stf dS( ~ I L dSi - Tm (2 - 2o) 

where T = 0, — 2n, or — 4w for points outside, on the boundary, or inside the volume 
defined by the surface S = Sb U So U Sl U Sr U Sp of Figure 2-1. 

The choice of Green function G is determined by the boundary conditions and the type 
of flow to be modeled. A complete development of the boundary integral equation method 
(BIEM) is available in Hunt [13]. 



The Combined High-order Spectral and Boundary Integral Equation Method 
The combined high-order spectral and boundary integral equation method (hereinafter re- 
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ferred to as the ’combined method’) eliminates the limitations of the base methods identified 
in Chapter 1. In particular, arbitrary body shapes are accommodated and all body bound- 
ary conditions are satisfied on the exact positions of moving bodies. The combined method 
defines a “total potential”, <E>x(x, y, 2 , t ), as the sum of a “spectral potential”, <J> 5p (:r, y, 2 , i), 
and a “body potential”, <J>b(z, y, 2 , t). These potentials are defined by the following rela- 
tionships: 



M 



$ T (x,y,z,t) = Y, 

m= 1 
M 

$s P (x,y,z,t) = Y 



m= 1 
M 



$B(x,y,z,t) = Y 



m=l 



( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 



In the combined method (z, y, 2 , £) is constructed to satisfy the initial-boundary value 
problem specified above. 

The function representing the body potential is chosen to set the body’s contribution to 
the potential on the mean free surface to zero. For a non-lifting body the potential is that of 
a periodic array of sources located on the body’s surface and its negative image. This choice 
is motivated by the sequence of boundary conditions used to determine the potential at each 
order, <J>( m )(x, 0, 2 , £), (equations 2.16 and 2.17). The sequence of boundary conditions to 
be satisfied on the mean free surface is 

(*» °> 0 = (* » °> 0 + (z> 0, 2, t) = $ s (x, 0, 2 , t) (2.24) 
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and 



(x, 0 ,z,t) = $ip : 1 (x, 0, 2 , t) + (x, 0, 2 , t) = 

m-1 i 

- £ Ti 7n * im - l) (x,0, z ,t) , m = 2, 3, . . . , M (2.25) 
^ /! dj/‘ 

As stated above, the body potential’s contribution to the velocity potential on the mean 
free surface is zero. With this choice the left hand sides of equations (2.24) and (2.25) are 
conditions to be satisfied by the spectral potential. At each order, the right hand sides 
are known Dirichlet boundary conditions. The coefficients of the Eigenfunction expansions 
which represent the spectral potential are determined directly from equations (2.24) and 
(2.25). 

In the combined method the body boundary condition is 



d$T 

dn 



d$sp 

dn dn 



— U ■ h on 5 b(Z) 



(2.26) 



where U is the body velocity and n is the body unit normal vector. 

In the combined method — d$ S p/dn\s B is added to the body’s normal velocity to yield 
the condition to be satisfied by the singularity distribution on the body. For a non-lifting 
body the condition is: 



d_ 

dn 




dj^?_ 

dn 



|s fl + U ' n ) 



m = 1 



(2.27) 



d_ 

dn 




dn 



|5b ) 



m = 2, 3, . . . , M (2.28) 
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where a(£) is the function describing the source distribution on the body surface. Equa- 
tions (2.27 and (2.28) can derived from the second form of Green’s theorem [17] or by 
directly assuming that the body can be represented by a distribution of sources along the 
body boundary. 

The coupling of the body and and spectral potentials occurs through the body bound- 
ary condition, equations (2.27) and (2.28), and the free surface boundary conditions, equa- 
tions (2.12), (2.13), (2.24) and (2.25) . The value of the body potential’s contribution to 
the total potential on the mean free surface is zero; the value of the vertical component of 
the gradient of the body’s potential on y = 0, d^B/dy\{x 1 o 1 z y t) ls n °t. The body contributes 
to the conditions to be satisfied on y = 0 through this quantity. 

Force 

The force on a body due to the dynamic pressure in an ideal, incompressible, irrotational 
fluid is obtained by integrating the pressure over the body surface [33]: 



F = J J pndS 


(2.29) 


M = [ [ p (r x n) dS 
J Js B 


(2.30) 



The useful forms of equations (2.29) and (2.30) are derived by applying to the unsteady 
form of Bernouilli’s equation either (1) the kinematic transport theorem and Gauss’ theorem 
or (2) the definition of the exact differential of the potential. If method (1) is applied the 
resulting force expression is [33] 



F - -4 IL * niS ~\ dS+ UL. (£) v * dS (2 - 31) 
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The following expression for the force is a result of the application of method (2) ([16] and 
[54]) 



F = -p- 7 - [ [ <j>ndS- l{V<t>- V0n) dS + p [[ ( U-V<fi)ndS (2.32) 

dt J J s q 2 J J s q 



The force expressions differ in the third term on their respective right hand sides. The 
difference between the expressions is explored in Chapter 5. 

Equation of Motion 

The equation of motion for a rigid body is derived from Newton’s Laws. The appropriate 
expression of dynamic equilibrium for a direct simulation of the free and forced motions of 
a submerged body is Newton’s Second Law [20] 

M X(f) = F(t) (2.33) 

M is the body mass matrix, X(t) is the body acceleration vector, and F(t) is the force 
vector. F(t) includes the wave and control surface forcing. The solution of the equation 
of motion for a body enables simulation of free and forced motions. The calculation of 
the wave forcing presently controls motion simulation computation speed and is a principal 
motivation for this thesis. 

A state variable representation of the body kinematics makes direct use of the solution 
of equation (2.33) to simulate body motions in a time varying flow 
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(2.34) 



/ t -tot « 

X(t) dt 

^ rt-\-6t 

<5X(f) = J X(t)dt 



(2.35) 



Equations (2.34) and (2.35) tire used to update the body velocity and position 



X(t + St) = X(t) + SX(t) 


(2.36) 


X(t + St) = X{t) + SX{t) 


(2.37) 



In the following section the equations used to implement the combined method in two 
dimensions are presented. 

2.4 Implementation of the Combined Method in Two Di- 
mensions 

The combined high-order spectral method - boundary integral equation method formulated 
for two and three dimensions has been implemented in two dimensions in this thesis. The 
purpose of this section is to present the equations used to carry out this implementation. 

The chosen representation of the body potential enables the initial condition for the 
free surface potential, = 0), to be used with equation 2.24 to determine the initial 

coefficients of the eigenfunctions which represent the spectral potential: 

$y } (x, 0, t = 0) = $£>(*, 0, t = 0) + $g ) (x, 0, t = 0) = $ s {x, t = 0) (2.38) 
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With the body potential represented as an array of sources and its negative image, 



4 1) (x,0,t)=0 

and the value of the spectral potential on the mean free surface is 



(2.39) 



$<#(x,0,t = 0) = $ s (x,t = 0) (2.40) 

The spectral potential eigenfunction expansion is represented as a Fourier series and the 
coefficients of each mode, $1^(0) axe determined by 

N 

*$(s.O ,t = 0) = £ *Si(0)’M*,0) (2.41) 

71=1 

The representation of the spectral potential is similar to that of the HOSM of Dommer- 
muth and Yue [8]. 

The initial body velocity, t/(0), and the gradient of the first order spectral potential, 
(x,y), are used with the first order body boundary condition to determine the first 
order source distribution on the body a^\x,y, t): 



8 _ 

dn 




x,y,t = 0)G(x, y, x\ y') <1 Sb 



ff - 5$sp(x,y) 

U - n ~ dn |SB 



(2.42) 



35 



where G{x, y, x y') is the potential of a periodic array of sources [34] and its negative image 
reflected about y = 0 



C(x, = ‘ log ( 2c osh;^-2cos^<^)) 

1. / , 7 t(x — x f> ) n n(y — (2h 0 — y')) \ . . 

- -log (2 corf, -2 cos ^ "» ) (2.43) 

where h 0 is the depth of the body’s centroid. 

The first order source distribution on the body and the first order spectral potential are 

the forcing for the solution of the second order problem through the boundary condition 
(2) 

imposed on 4>^(:r, 0, £): 



4 2) (x,0,t) = $$(X,<M) + * { £ ) (x,0,t) (2.44) 

( 2 ) 

As with the first order body potential, (x,0,£) = 0- Therefore equations (2.44) and 

( 2 ) 

equations (2.25) are combined to yield the condition to be satisfied by $ K S p(x,0,t): 



*<2) (I o 0 - _ asjj'fro ■*) 



s(D/ 



(2.45) 



The second order spectral potential is used in the body boundary value problem to solve 
for the second order body potential. 



dn 






d$% } 

dn 



Isb 



(2.46) 



The sequence used to determine the second order source strengths and the higher order 
spectral potentials follows that used for the first order quantities above. 
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With the body and spectral potentials known to the desired order M, the quantities 
needed to compute the pressure and force on the body are known. For a body free to 
respond to the waves, the force is computed and used to update the body’s state using the 
methods described above. 

The body and spectral potentials are used with the Zakharov equations to update the 
free surface potential and free surface shape. The quantities needed are 



dy{x,t) d$ s (x,t) d$ T {x,ri,t) /0 

ST’ — '“ d Ty (247) 

dr)(x,t)/dx and d$ s (x,t)/dx are calculated in spectral space, by the method used by 
Dommermuth, Yue, and Liu [23]. d<&T{x,7],t)/dy IS calculated by determining the contribu- 
tions of the body and spectral potentials separately. The spectral potential’s contribution 
is computed as in [23] and [8]: 



d* ap (x,T,,t) d k+1 



dy 



E E frE< , < t >3 ? +T > M*,0) (2.48) 

7i—l k=0 * n=l * 



The body’s contribution at each order to the fluid vertical velocity on the free surface, 



dy 

can be computed directly on the free surface or by Taylor series expansion of the vertical 
component of partial derivatives of the body potential evaluated on y = 0. For k > 1, 
d k ^/dy k is made computationally efficient using Laplace’s equation to reduce the number 
of y derivatives required [23]. 
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Chapter 3 



Numerical Implementation 



The initial boundary value problem defined in the previous chapter is solved numerically. 
The purpose of this chapter is to describe the numerical implementation of the mathematical 
formulation. This chapter will: 

• Present the problem in the form of systems of discretized equations to be solved by 
appropriate numerical methods; 

• Analyze the accuracy, consistency, and stability of the proposed numerical scheme; 
and 

• Demonstrate necessary conditions for the method to converge to solutions of the 
physical problems of interest. 

3.1 System of Equations 

I 

The method developed for solving the initial boundary value problem requires the solution 
of a system of equations at each time step. The system of equations is determined by the 
specified form of the boundary conditions in the initial boundary value problem. The form 
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of the free surface boundary conditions, the periodic nature of the problem formulation, and 
the choice of a boundary integral representation of the body determine the characteristics 
of the numerical formulation. 

The numerical solution is enforced at points along the boundary of the computational 
domain S = S^U So U SlU SrU Sp of Figure 2-1. The solution is enforced at Nx evenly 
spaced points along the horizontal (x) axis, Nq evenly spaced points around the body, and 
No evenly spaced points along the bottom. 

The system of equations developed to represent the boundary value problem at each 
time step is expressed symbolically in matrix form as 



[F] [Lbw] [L„„] 

[Fb] 1 L bb] [Lob] 

K] [LboI [Lool 





’ [*!?'] 






«4 m) 


== 




*i> m) 








. 



[R( m >] 



U • n | Ss for m = 1 



0 for ra>2 

[0] 



F is the array of eigenfunctions for the spectral potential: 

e ikoX\ e \ko\yi gikixi e \ki\yi e ifc N x Xl e^ N x^ yi 

e ikoX2 e \ko\y2 e ik\X2 e \k\\y2 ... e^N x X2 ^N x \ y2 

[F] = 

}koX Nx |fco|y Nl lk\X Nx |fcl|y Nl . .. e lk N z X N x e \ k N z \ y N 



(3.1) 



(3.2) 
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F' B ] represents the influence of the spectral potential on the body: 



ft] = (v( ) • «is. 



(3.3) 



• F 0 represents the influence of the spectral potential on the bottom: 



ft] = m i • fi] Sc 



(3.4) 



[L bw ] defines the influence of the body potential on the spectral potential on y=0 



[L ow ] defines the influence of the bottom potential on the spectral potential on y=0 



[L bb ] defines the body’s self-influence 



[L ob ] defines the influence of the bottom potential on the body 



[L bo ] defines the influence of the body potential on the bottom 



[L 00 ] defines the bottom’s self-influence 



$sp (m) ] is the array of spectral potential modal amplitudes, A n (t) 



is the body potential 



$ 0 ^ is the bottom potential 



Note that [L 00 ] and [L ow ] need to be computed only once. 
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[R (m )J is the function representing the Dirichlet boundary conditions for 



rU) = $5 

R<m) = (3.5) 

where (x, 0, t) = (f)^~ k ^ (x, 0, t) + <j)^~ k ^ (x, 0, t ) + <f^o~ k ^ (x, 0, t). 

Note that all the [L ] f s contain the influences of the sources and their negative images. 

For deep water, 3.1 is rewritten as: 



[F] 


l^BW 


K1 


it-..: 







for m = 1 



(3.6) 



L JL J 0 for m > 2 

The body potential is represented as a horizontally periodic array of sources with neg- 
ative images reflected about the mean free surface, y = 0. When present, the bottom 
can be similarly represented. A consequence of this choice is that the [L BW ] and [L ow ] 
matrices are identically zero. It is noted that this construction may appear similar to the 
’pressure release’ problem, the high-frequency limit of the linear free surface boundary con- 
dition (cf. Newman [33]). However, in that linearized problem the complete solution yields 
4>(x, 0, t) = 0 whereas this is not the case for the problem at hand. The negative images 
here are used to provide computational efficiency without restricting the generality of the 
solution. In particular it permits fast calculation for the spectral potential at each time 
step. 
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3.2 Computational Steps 



The time simulation of wave-body interactions, developed here for the deep water case for 
simplicity, proceeds as follows: 

1. Initial conditions ( t = 0) are specified for: 

The free surface profile, 77 ( 0 :, ^); 

Free surface potential, $ s (a;,£); and 

Body position and velocity, and U(x,y,£),. 

2. The boundary value problem is solved for the spectral potential, (f>sp \ and the source 
strength distribution on the body, . For the first set of computations, at m= 1: 

a. The spectral potential is determined by solving the following for 



[F] [*£>] = [R (m) ] 



(3.7) 



b. With the spectral potential at a given order known, the source strengths are 
determined by solving 



[ F b] [*£>]+ [I'm] [*B m) ] 



0-S\s B 


for 


m = 1 


0 


for 


m > 2 



(3.8) 



The above steps (2. a and 2.b) are repeated up through perturbation order M. The 
body contributes to <j > through its vertical derivatives. 
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3. Compute the quantities necessary for the time integration of the free surface boundary 



conditions. These quantities are: 
and 

Vx (x,t). 

4. Integrate the free surface boundary conditions, equations (2.12) and (2.13), to update 
the free surface potential, 4> s (:r, t), and free surface elevation, 1 i(x,t) 

Steps 2, 3, and 4 are repeated at each time step. 

The following section details the manner in which each of the above steps is performed. 

3.3 Solution of the Boundary Value Problem 

Initial Conditions- Free Surface Elevation and Potential 

The solution of the nonlinear wave body interaction problem requires the free surface 
elevation and potential to be accurately specified. The free surface profile and velocity 
potential used here for initial conditions are generated by Stokes’ expansions. The mapped 
equations of Schwartz [41] are used to specify a Stokes’ Wave of arbitrary order N [22]. For 
the typical wave-body interaction analyzed in this thesis N = 20 is sufficient for specifying 
the initial free surface profile and velocity potential. The differences in the free surface 
profile for expansions of order 2, 10 ,and 20 are shown in Figure 3-1. 
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Free Surface Elevation 



Free Surface Profiles by Stokes’ Expansion 




Figure 3-1: Initial free surface profile 
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Though the calculations done here are for Stokes’ waves as initial conditions, the equa- 



tions and solution methods apply to fully arbitrary initial conditions. 

Influence Coefficient Matrices: [F] and |V B j 

The spectral potentials at each order, are obtained by solving equation 3.7 for the 
Fourier modal amplitudes, A n (t). In expanded form: 



N x 



ct>W(xAt) = Y, A n(tW k '' X = 

n = 0 

$ s (x,t) 



m = 1 



k - ElfcLi 1 [$■ ($p k) (*» y » 0 + 0b* k) (*. y . <)) ] y=Q m > 2 



(3.9) 



The higher order spectral potentials are functions of the free surface elevation and 
vertical derivatives of the lower order potentials evaluated on the mean free surface. The 
vertical derivatives are obtained by computing the modal amplitudes of the lower order 
potentials by Fast Fourier Transform (FFT), multiplying the modal amplitudes by |fc n |, 
and generating the derivatives in physical space by inverse FFT: 



” 71 = 1 



(3.10) 



The products needed in the boundary conditions for the spectral potentials, e.g., 



f?(z)^-(z,(M) 



(3.11) 



are formed in physical space. The operational count associated with obtaining each is 
0{N X log N x ). 
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With the Fourier amplitudes of the spectral potentials known, the spectral potential and 
its derivatives can be determined throughout the fluid domain. In particular, the influence 
of the spectral potential on the body is computed directly from 




A 0 (t) [(ike ik ° x ')i + (Ne' fc °l*)j) • n B 
Ai{t) ((ike iklX ')i + • n B 

A 2 {t) ((ike lkiXi ) i + (Ifolel* 0 ^)]^ • n B 



(3.12) 



A Nx ((ike ik »* Xi )i + 






Xj,yi are the locations of the Nb body segment midpoints. 

The operational count associated with determining the influence of the spectral potential 
on the body is 0(N x Nq )• The influence of the spectral potential on the body is combined 
with the body boundary condition to form the right hand side for the equation to be solved 
for the source distribution on the body: 



[a(xi,y u t)] = [L bb ] 1 




(3.13) 



Equation (3.13) is solved by factoring the body-body influence coefficient matrix, [L BB ] _1 , 
and determining the source strengths by back substitution. The body-body influence co- 
efficient matrix must be factored only once, at the first time step. The operational count 
associated with solving 3.13 at the first time step is O(Ng). The operational count at 
subsequent time steps is O(NbNb). 
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Equation 3.13 requires derivatives of the form 



a fc «4 m - fc) 

d k y 



(i,0, t), 



k = 1, 2, . . . , m — 1 



(3.14) 



50^(x,O, t)/dy is computed using the analytical expression for the representation of the 
body as a distribution of periodic sources and their negative images. Higher derivatives 
are found through the use of Laplace’s equation ([8]). For example, /<9y 3 ) | y=0 IS 

required for the solution of the fourth order problem. d(f)^ /dy is a harmonic function, 
therefore 



= (^£\ 

dy 2 \ dy ) dx 2 \ dy J 



(3.15) 



/dy is computed as before and the horizontal derivatives are computed numerically. 
The symmetry properties of the body-image system force the even numbered vertical deriva- 
tives to be zero. 

After the initial set up effort, the total operational count associated with solving the 
boundary value problem to order M at each time step is 



0{M[N x log (iV x ) + N x N b + N b N b }) (3.16) 

The total operational count of the High Order Spectral Method (HOSM) of Liu, Dom- 
mermuth, and Yue [23] is 



0(M[N x log (N x ) + N x N b + N b log N B }) (3.17) 



47 



As in the HOSM, typically N x 3> Nb- Thus, for large numbers of panels the new method 
is slower than the full spectral method, but provides the benefit of allowing completely 
arbitrary body geometries. 

3.4 Updating the Free Surface Potential and Elevation 

The Zakharov equations (equations 2.12 and 2.13) are used as time evolution equations 
for the free surface elevation r)(x,t ), and the free surface potential The time 

integration requires (a) the accurate determination of the spatial derivatives of the free 
surface elevation r ] x , free surface potential and velocity potential </> y , and (b) a stable 
and accurate time integration scheme. 

Horizontal Derivatives of r](x,t) and <f) S (x,t) 

The free surface elevation and free surface potential are defined at evenly spaced points 
along the x axis. The horizontal derivatives of these quantities are most efficiently performed 
in Fourier space: 



= £ B n (t)ik n e iknX (3.18) 

71 = 0 

(3.19) 

71=0 

where B n (t) and A^(t) are the Fourier amplitudes of the free surface elevation and potential 
respectively. The operational count associated with computing these spatial derivatives is 
0(N X log N x ). Products involving these quantities are computed in physical space. 

The Vertical Velocity of the Fluid on the Free Surface 
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The vertical velocity of the fluid on the free surface is the sum of the contribution from 
the spectral potential and the body potential. The contribution from the spectral potential 
is determined by computing the derivatives in Fourier space [8]: 

^fc+T« r »(*,y)ly=o = IM fc+1 *»(*>y)U (3.20) 

The Fourier amplitudes are found by FFT and the appropriate power of the wave number 
is multiplied by the series representation of the function, and the required derivative of the 
function is computed by the inverse FFT. The operational count associated with computing 
these derivatives is 0(N X log N x ). 

The contribution, at each order, of the body potential to the fluid vertical velocity on 
the free surface is computed using the derivative of analytical expression for the potential 
due to the body and its negative image. 

The operational count associated with computing the vertical derivatives is 0(N x Nb)- 
The total operational count associated with computing the terms required for solution of 
the kinematic and dynamic free surface boundary conditions is 0(M[N x \ogN x + N x Nq + 
N B N E ,]). 

The above steps provide the information required to compute dr)(x, t)/dt and d<fi s (x, t)/dt 
at each time step. In the next section the method used to integrate the free surface evolution 
equations in time is presented. 



Time Integration 

The integration scheme used to evolve the free surface shape and free surface potential 
is a fourth order Runge-Kutta (RK4) scheme. The stability properties of this method as 
applied to the non-linear wave body interaction problems were developed by Dommermuth 
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[7] and used by Liu [22] in his extension of Dommermuth’s work. The important points 
from Dommermuth are: 

1. RK4 is conditionally stable; 

2. The Courant condition for RK4, based on von Neumann stability analysis of the 
linearized free surface conditions, is g k gr i(iAt 2 < 8, k gr id = 2ir/Ax ; and 

3. RK4 is mildly dissipative. 

Analyses of other time integration schemes for numerical simulations of wave-wave and 
wave-body interactions are found in Kring [20] and Dommermuth [7]. 

Considerable information is available concerning the development and growth of “saw- 
tooth instabilities” in time simulations of free surface flows. Olmez [37] and Dommermuth 
[7] provide comprehensive background discussions of the saw-tooth instability. Longuet- 
Higgins and Cokelet [25] were the first to report the presence of the sawtooth instability 
and contended that its presence was partly physical. Other investigators reporting sawtooth 
type instabilities are Baker et al. [1], Lin et al. [21], and Sen et al. [42]. 

Dommermuth [7] and Roberts [39] stated that the sawtooth instabilities are numerical. 
Olmez [37] stated that the cause of the instability is unknown. 

The discretization of the physically continuous problem introduces a sawtooth wave 
form at the highest resolvable wave number, k = 2n/ Ax (Figure 3-2). 

The coupling provided by the method of solution outlined in the previous chapter ensures 
that all modes at all orders are coupled. Instabilities at any mode are coupled to all other 
modes. The presence of physical or numerical stabilizing mechanisms determines whether a 
given frequency is stable. In the chosen mathematical formulation of the physical problem 
the fluid is assumed to be ideal. Therefore no physical dissipative mechanism exists. The 



50 



The value of the eigenfunction at the highest wave number 
e (lkMA * > , k=rc/Ax 




Figure 3-2: Sawtooth at Highest Wavenumber 

numerical scheme must provide the dissipation absent in the chosen physical model. RK4 
is mildly dissipative. 

The level of accuracy of the spatial derivatives evaluated on the free surface influences 
the onset of sawtooth instabilities (Roberts [39] and Dommermuth [7]). Derivatives of the 
potential of each order m are of the form 

= (3.2i) 

where n is the Fourier mode. The magnitude of the derivative is proportional to ( k n ) m . The 
error at a given mode and order is amplified by this factor [7]. If dissipation is not provided 
the growth of this error can be unbounded. Longuet-Higgins [25], Dommermuth [7], and Liu 
[22] all use a form of smoothing to further control high wave number instabilities. In this 
thesis an ideal low pass filter is used to remove the highest wave number modes. The cutoff 



51 



wave number is constrained by physical considerations. The filter must allow physically 
relevant waves to develop, A C utof f / Dbody 1 and must not remove appreciable energy from 
the system. The filter is implemented by transforming the free surface elevation and free 
surface potential into the frequency domain by the FFT. All modes above the cutoff wave 
number are removed, and the filtered free surface elevation and potential are generated by 
inverse FFT. 



3.5 Convergence 

The Lax-Richtmyer equivalence theorem states that a consistent approximation to a well- 
posed linear problem is stable if and only if it is convergent [10]. First, it is necessary to 
demonstrate that the numerical solution method is consistent. The consistency is demon- 
strated by conservation of volume and conservation of energy. Confidence in converged 
results is provided by comparison of converged results with experimental results and results 
generated by other numerical methods (Cointe [6]) . 
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Energy and Volume Conservation 



Volume conservation requires 




dx = Constant 



Energy conservation requires 



»JJL G** +»)<«'-* 



(3.22) 



(3.23) 



A convenient form of the total energy can obtained by making use of the kinematic infor- 
mation, U = V(f> and V 2 ^> = 0 and the first form of Green’s theorem 



[(j>V 2 <j> + (V<£) 2 ] dv = / n • (j)V(j)ds 



(3.24) 



The above information can be used to express the kinetic energy T as (Milder [28] and 
Kinsman [18] 



T = \ p J v< ^ ■ V<t)dv = \ p / w ' ds = \ p I ^ Srit dx (3.25) 



The expression used for computing the total energy £, in the system is 



C = \p / dx+^pg J T) 2 dx 



(3.26) 



The total energy is computed at each time step. 
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Convergence- Discussion 



The convergence properties of the combined HOSM-BIEM method are functions of the 
following sources of error [7] : 

1. Finite numbers of body segments Nb, Fourier modes N x , and perturbation order M ; 

2. Time integration scheme; 

3. Finite computational domain; 

4. Filter cutoff wave number Nfu ter ; and 

5. Aliasing. 

The first four sources of error require no explanation. The convergence properties associated 
with each will be examined in the subsequent section. 

Solutions associated with wave numbers above the highest grid wave number ( k = 
27r/Ax) can not be distinguished from waves of lower frequency [46]. This is aliasing and 
it occurs is when products are formed in the spectral method, e.g., r] x rj x 

£ %(,)«*■*(*» £ %(,)«“■*(*» = -kl fa.W) 2 e 2i *-' + ■ • • + (»„,) 2 e“«. (3.27) 

n=0 n — 0 

Products are made alias free using the method developed by Orszag [38]. 
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Convergence- Demonstration 



The convergence properties of the HOSM-BIEM are demonstrated in this section. 

Nomenclature: 

• M = Perturbation order 

• Nb = Number of body segments 

• F y = Steady vertical force 

• p = Fluid density 

• g = Gravitational constant 

• A = Wave amplitude 

• k = Wave number 

• R = Body radius 

• H = Distance from mean free surface to center of body 

• Nw = Domain size; the number of waves of the fundamental wave length in the 
domain 

• N x = Number of Fourier modes; equal to the number of segments along the x axis 

• T = Period of the fundamental wave 

• Ts = Simulation time 

• A t = Time step size 

• N filter = Fourier mode for low pass filter cut-off 
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N b 


M = 2 


M = 3 


M = 4 


64 


0.2807 


0.2799 


0.2798 


128 


0.2783 


0.2775 


0.2775 


256 


0.2768 


0.2760 


0.2760 



Table 3.1: Convergence of the vertical drift force, F y /pgA 2 , with body segments Nb and 
order M . kA — 0.04, kR — 0.4, and H/R = 2; and N\v — 16, N x — 64A/\y, iV filter = 24iVvv 
for N b = 64 and N B = 128, Njuter = 32 N w for N B = 256, T/At = 128, T 5 = 5T. 

Table 3.1 shows the convergence of the vertical drift force with increasing number of 
body segments Nb and perturbation order M. The convergence rate appears to be better 
than linear with increasing number of body segments, and exponential with increasing 
perturbation order. The full spectral method of Dommermuth converges exponentially 
with increasing number of body modes (segments). The combined method of this thesis 
uses a boundary integral representation for the body’s contribution to the velocity potential. 
Therefore the convergence rate is not expected to match that of the full spectral method. 
The decreased convergence rate is a penalty associated with the improved flexibility of the 
present method. 

Figure 3-3 shows the behavior of the solution for the vertical force as the value of the 
low-pass filter cut-off wave number is varied. The cut-off wave number corresponds to waves 
whose lengths vary from 2/3 to 1/48 of the wavelength of the fundamental of the initial 
Stokes wave. For cut-off wave numbers corresponding to waves shorter than 1/4 of the 
fundamental, there is no appreciable change in the vertical drift force. 

Table 3.2 shows the convergence of the vertical drift force with increasing number of 
Fourier modes and perturbation order. The convergence rate appears to be exponential with 
increasing number of Fourier modes and exponential with increasing perturbation order. 
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0.30 r 




Figure 3-3: Behavior of the vertical drift force, F y /pgA 2 , with varying cut-off wave number 
N filter • kA = 0.04, kR = 0.4, and H/R = 2; and Nb = 256, Nw = 16, N x = 64Nw, 
T/At = 128, M = 4, T s = 3T. 



Table 3.3 shows how the convergence of the vertical drift force varies with increasing 
wave slope. The convergence rate with increasing number of Fourier modes decreases as 
the wave slope increases. N x = 32 Nw Fourier modes was insufficient for kA = 0.12. The 
number of Fourier modes was insufficient to resolve the local steep waves in the vicinity of 
the body and the numerical solution broke down after 2 wave periods. 

Table 3.4 shows the convergence of the method with respect to time step size. The 
behavior of the solution with time step size is consistent with the findings of Korsmeyer 
[19]. For solutions that are converged with respect to the spatial variables, there is a wide 
range of solutions for which changes in time step size do not affect the accuracy of the 
solution. For larger time step sizes the solution error grows rapidly and is unstable [19]. 
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€ 


N x 


M = 2 


M = 3 


M = 4 


0.04 


32 N\y 


0.2769 


0.2758 


0.2760 


64 N w 


0.2768 


0.2760 


0.2760 


128 Nw 


0.2768 


0.2760 


0.2760 



Table 3.2: Convergence of the vertical drift force, F y /pgA 2 , with Fourier modes N x and 
order M. kR = 0.4, and H/R = 2; and Nw = 16, Njut er = 24Nw> Ng — 256, T/ At 128, 

T S = 5T. 



€ 


32 N w 


64 N w 


128 Nw 


0.04 

0.08 

0.12 


0.2760 


0.2760 


0.2760 


0.2760 


0.2706 


0.2706 


- 


0.2692 


0.2640 



Table 3.3: Convergence of the vertical drift force, Fy/pgA?, with Fourier modes N x as wave 
slope kA increases. M = 4, kR = 0.4, and H/R = 2; and Nw = 16, N/u ter = 2ANw, 
N b - 192, T/At = 128, T s = 3 T. 



Figure 3-4 shows the behavior of the HOSM-BIEM method as the domain size is varied. 
The purpose of evaluating the behavior of the solution as the domain size varies is to 
determine the influence of the periodic problem formulation on the ability to achieve and 
sustain limit-state behavior before “reflections” from the periodic boundaries effect the 
solution. Figure 3-4 shows that limit state behavior is reached within 3 periods, Ts = 3T, 
of simulation time. The effect of “reflections” for Nw = 8 is evident in Figure 3-4. 
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T/At 


II 


M = 3 


M = 4 


32 


0.27826 


0.2775 


0.2774 


64 


0.2783 


0.2775 


0.2775 


128 


0.2783 


0.2775 


0.2775 



Table 3.4: Convergence of the vertical drift force, F y /pgA 2 , with time step size T/At and 
order M. kA = 0.04, kR = 0.4, and H/R = 2; and Nw — 16, Nfuter = 24 Nw, Np — 64Nwi 
N b - 128, T s = 5T. 




Figure 3-4: Behavior of the vertical drift force, F y /pgA 2 , with varying domain size Nw 
where Nw is the ratio of the domain size to the fundamental wavelength. kA = 0.04, 
kR = 0.4, and H/R = 2; and Nb = 128, N x — 64Nw, N/uter = 24 Nw, T/At = 128, 
M — 4 
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Chapter 4 



Wave Diffraction by Submerged 
Cylinders 



The purpose of this chapter is to present the results for the analysis of the diffraction of 
incident waves by submerged circular cylinders. The diffraction of waves by cylinders of 
circular and ’pontoon’ cross sections are studied. The pontoon case is analyzed to demon- 
strate the applicability of the method to arbitrary geometric shapes. Figure 4-1 defines the 
problem parameters. 

The oscillatory and steady forces for a range of problem parameters are presented and 
compared with related results generated from the other solution methods and experimental 
results discussed in Chapter 1. 



4.1 Diffraction by a Cylinder of Circular Cross Section 

The forces resulting from the interaction of a fully submerged circular cylinder with Stokes’ 
waves are analyzed in this section. The forces on the cylinder are evaluated as functions of: 
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Figure 4-1: Diffraction problem parameters 

• Keulegan-Carpenter number, Kc = ne~ kH A/R ; 

• Wave slope, kA ; and 

• Submergence ratio, h = H/R. 

Although a “natural” nondimensional variable is kR, the variable Kc, which is equiv- 
alent in the presence of the other nondimensional variables, is used instead to make the 
results directly comparable to others in the literature. 

The forces analyzed are: 

• First, second, and third harmonics of the vertical and horizontal forces, and 

• Steady vertical and horizontal forces. 
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Method of Computing Forces 



Time simulations of the initial boundary value problems are performed long enough 
for the solution to achieve limit-state behavior and terminated when “reflections” from the 
domain boundaries begin to influence the solution. Limit state behavior is typically reached 
within two periods. The length of time the limit-state behavior is sustained is determined 
by domain width, Nw (see figure 3-4). For the diffraction problem the forces are computed 
by applying Equations 2.29 and 2.30 to the time history of the sum of the potentials at 
each order, , m = 1,2,..., M, to generate a force time history. For stationary bodies 
the equations reduce to 



F{t) = -p[ p(t)ndS = -p[ 

Js B Js b L 






ndS 



(4.1) 



where 4>r(t) is the total potential on the body’s surface. 
For the circle Equation 4.1 becomes 



F(t) = -p [ 

Js B 



dt 2 \r 0 dd ) 



ndS 



(4.2) 



where r D is the circle radius. 
Numerically, 



A r B 

m = -p'Z 



i=l L 



d<f>i T (t) 1/1 d<t>i T {t) 



dt 



+ 



1 /i d<f> iT (t) v 

2 \r 0 de ) 



n idS 



(4.3) 



d<pi T (t)/dt and (d<j) T (t)/dd) i are computed by temporal and spatial central differencing, 
respectively. The potential and pressure along each segment are assumed to be constant. 
Therefore the integrated quantity at each time step is computed as the simple sum, 
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(4.4) 



N b 

i = 1 

The force at each order m is related to the n th harmonic coefficient of the force: 

F n = F (m) ~ i [ t+T F{t)e> mult dt-, n — m — 1,2, ...,M (4.5) 

T Jt 

with j = \f—\- 

The forces computed by the combined method are compared with solutions generated 
by the other methods described in Chapter 1. As noted by Liu, et al. [23], there is no direct 
relationship between the terms of the nonlinear solution and the F n = e m 4>^ 

terms of a perturbation series solution generated by a boundary condition perturbation 
frequency-domain approach. Liu et al. demonstrate that relationships between the nth - 
harmonic of the nonlinear solution, = J2m=i n = 1, 2, ..., N x are related to the 

frequency domain amplitudes, by 




i? = K 2> + 0(< 3 ) 


(4.6) 


1 = 4( m> + 0(«" +1 ) 


(4.7) 



Numerical Results 

In Figure 4-2 the first harmonic of the horizontal diffraction force on the submerged 
circular cylinder is shown as a function of the Keulegan-Carpenter number. The results are 
presented together with: 

• The experiments of Chaplin [3] and Miyata [30]; 
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• The linear theory of Ogilvie [35]; 



• The High-order Spectral Results (HOSM) of Liu, Dommermuth, and Yue [23]; and 

• The Mixed Eulerian-Lagrangian (MEL) method of Cointe [6]. 




Figure 4-2: The first-harmonic of the horizontal diffraction force. Experiment (x); linear 
result( — ); HOSM(CD); MEL(A); present method (•). ( kR = 0.21, H/R = 2.0.) 



The agreement between the method developed in this thesis and the other numerical 
methods for this case helps to confirm the theoretical and numerical validity of the new 
method. 

In Figure 4-2 the experimental results differ appreciably from the theoretical and nu- 
merical results at values of the Keulegan-Carpenter number greater than 0.5. The difference 
is due to flow separation, a fluid behavior not accounted for in the potential flow equations 
used in each of the theoretical and numerical methods cited. 
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The lack of agreement between the experimentally generated results and the results 
generated by the numerical methods is due to lift oscillating at the fundamental frequency 
in antiphase with the inertia force (Chaplin [3]). Chaplin’s results were consistent with the 
analytical results of Longuet-Higgins [24] who computed the magnitude of the circulation 
associated with the oscillatory flow around the cylinder. The decrease in vertical force 
with increasing Keulegan-Carpenter number is due to the (negative) lift asscociated with 
the circulation. 

The zeroeth, second, and third harmonics of the horizontal force are shown in Fig- 
ures 4-3, 4-4, and 4-5 along with the results generated by experiments, second-order theory 
(Ogilvie [35]), perturbation theory (Vada [48]), and HOSM and MEL computations. The 
results of the theoretical and numerical methods agree with the experimental results over 
a wide range of Keulegan-Carpenter numbers. 




Figure 4-3: The steady vertical diffraction force. Experiment (x); MEL(D); second-order 
theory( — ); present method(*). (kR = 0.21, H/R = 2.0.) 



65 




Figure 4-4: The second-harmonic of the horizontal diffraction force. Experiment (x); 
HOSM(D); MEL(A); perturbation theory( — ), present method(«). ( kR = 0.21, H/R = 
2 . 0 .) 
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Figure 4-5: The third-harmonic of the horizontal diffraction force. Experiment (x); 
HOSM(D); present method(*). ( kR = 0.21, H/R = 2.0.) 
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The drift forces as functions of wave slope are shown in Figure 4-6. The results of the 
present method are compared with those generated by the HOSM for the vertical force 
and the HOSM and MEL for the horizontal force. The normalized vertical drift force is 
independent of the wave slope. The normalized horizontal drift force is shown to decrease 
in magnitude as wave slope increases. The negative horizontal drift force was studied 
extensively by Dommermuth [7] and Liu et al. [23] and is not a principal focus of this 
thesis. Cointe [6] also computed a negative horizontal drift force. In Cointe’s work the 
negative horizontal drift force is only present for wave steepnesses close to the limiting 
steepness of the method. 

As wave slope becomes vanishingly small, the fourth order force should approach a 
constant. As seen in Figure 4-6 none of the numerical methods capture this result. The 
errors, however, for the high order spectral method and the present method are small. 




Figure 4-6: Drift forces as functions of wave slope, kA. F y is the vertical force; F x is the 
horizontal force. HOSM(D); MEL(A); present method(x). ( kR = 0.4, H/R = 2.0.) 
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The dependence of the drift forces on submergence, H/R, is shown in Figures 4-7 and 4- 
8. The experimental results are from Miyata [30]. As the body approaches the free surface 
the effects of wave breaking become important and the present method can not account 
for overturning waves. This limitation also exists with the MEL method of Cointe and the 
HOSM of Liu et al. 




H/R 



Figure 4-7: Vertical drift force as a function of submergence H/R. Experiment (x); 

HOSM(D); linear result( ); present method(O). ( kR = 0.4, kA = .12 for HOSM 

and Experiment; kR = 0.4, kA = .08 for present method) 
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Figure 4-8: Horizontal drift force as a function of submergence H/R. Experiment (x); 
HOSM(D); present method(O). ( kR = 0.4, kA = .12 for HOSM and Experiment; kR = 0.4, 
kA = .08 for present method) 
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4.2 Diffraction by a Cylinder of Pontoon Cross Section 



The ability of the combined method to compute forces on bodies of arbitrary geometry 
is demonstrated by computing the forces resulting from the interaction of a cylinder of 
pontoon cross section with incident waves. The pontoon shape is shown in Figure 4-1. The 
cylinder’s geometric definition is the one used by Vada [48]: 



x — 



y=b 



cos Q — a cos 3 Q 
1-a 
sin 0 + a sin 3 9 



l—o 



For the shape evaluated in this thesis, o = 0.1 and b = 0.5. 

Method of Computing Forces 



For the pontoon Equation 4.1 becomes 



(4.8) 



F(t) = -p [ 

Js, 



[ 


d<P T {t) 1 




ls B 


dt 2 1 


\ ds J _ 



ndS 



(4.9) 



Numerically, 



N b 

F(t) = -p'£ 

i = 1 

d(pi T (t)/dt and ( d(f> T (t)/ds) i are computed by temporal and spatial central differencing, 
respectively. 

The harmonics of the force time histories are determined using the method presented in 
the previous section. 



d<f >i T (t) 1 f d<M*) V 

dt 2 V ds ) 



n; dS{ 



(4.10) 
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Numerical Results 



Time simulations of the interaction of the pontoon with waves were performed and the 
resultant forces compared with those of Vada[47]. The zeroeth, first, and second harmonics 
of the horizontal and vertical forces as functions of the body submergence, h — H/R, and 
nondimensional wave number, K = R/(oj 2 g) are presented in Figures 4-9, 4-10, and 4-11. 
The results for the two methods are in agreement for the range of depths and wave numbers 
evaluated. 




0 2.0 40 

K 



Figure 4-9: First- and second-order horizontal oscillatory forces on the pontoon, (a) H/R = 

1.0; (b) H/R = 1.25; (c) H/R = 1.5; Vada ( , \F x i/pgRA\; — , \F X 2 /pgA 2 \)\ present 

method (D,\F xi /pgRA\ ; 0,\F x2 /pgA 2 \) 
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Figure 4-10: First- and second-order vertical oscillatory forces on the pontoon, (a) H/R = 

1.0; (b) H/R = 1.25; (c) H/R = 1.5; Vada ( , \F y \/ pgRA\‘, — , \F y 2 / pgA 2 \)\ present 

method (□, \F yX /pgRA\ ; O, \F y2 /pgA 2 \) 
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I 



Figure 4-11: Steady vertical force on the pontoon, (a) H/R = 1.0; (b) H/R = 1.25; (c) 
H/R = 1.5; Vada (— , \F y / pgA 2 \)-, present method (□, \F y \/ pgA 2 \) 
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4.3 Summary 



The ability of the combined high-order spectral- boundary integral equation method to 
estimate the forces due to the interaction of a stationary two dimensional cylinder of arbi- 
trary cross section with waves has been demonstrated. The method enables second order 
and higher forces to be computed with accuracy consistent with the accuracy of high-order 
spectral and Mixed Eulerian-Lagrangian methods for similar problem parameters. The 
present method computes free surface quantities faster than an MEL and can accommodate 
more general body shapes than the high-order spectral method. 

In the next chapter, the combined method will be applied to the radiation problem for 
a fully submerged cylinder undergoing large amplitude forced oscillatory motion. 
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Chapter 5 



Wave Radiation by a Circular 
Cylinder Undergoing Large 
Amplitude Oscillatory Motion 



The purpose of this chapter is to demonstrate the application of the combined high order 
spectral-boundary integral equation method to the problem of a fully submerged circular 
cylinder undergoing large amplitude forced oscillatory motion beneath an initially quiescent 
free surface. Figure 5-1 defines the problem parameters. 

The cylinder is forced to oscillate in heave and in circular orbital motion. The free 
surface profiles and the oscillatory and steady forces resulting from the forced motions are 
presented and compared with results generated by the other solution methods described in 
Chapter 1. 
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Figure 5-1: Radiation Problem Parameters 

5.1 Radiation Problem for a Circular Cylinder in Heave 

The forces and free surface profiles resulting from the interaction of a fully submerged 
cylinder forced to oscillate in heave with an initially quiescent free surface are presented in 
this section. The forces are evaluated as functions of: 

• The nondimensional oscillation amplitude , £/i?; 

• The submergence ratio, h = H/R\ and, 

• The body radius to wavelength ratio (x 27 t), kR = 2irR/\. 

The forces analyzed are: 

• The first and second harmonics of the vertical force; and, 

• The steady vertical force. 
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Method of Computing Forces 



The mathematical formulation developed in this thesis solves the initial boundary value 
problem in an earth fixed inertial reference system. For the radiation problem the expression 
for the force is developed from Equation 4.1: 



FW ->IL 



d<t> T {t) 1 



dt 






ndS 



(5.1) 



Where (pr is the total potential on the body surface. 

The chosen numerical formulation provides information to compute d(pr/dt\ SB directly. 
However, the body is moving so that dcpr/dt ^ d<pr/dt. The exact differential of the 
potential is used to develop an expression for d<pr/dt: 



d<pT d(pT dx d(f>T dy d<pr 

dt dt dt dx dt dy 



d(pT 

dt 



— U • V(pT 



(5.2) 



Equation (5.2) is combined with Equation (5.1) to provide the expression used for com- 
puting the force on the moving body [16]: 



F = 




- (V<£r • V<£r) + 



(U-V<fr) 



n dS 



(5.3) 



In Chapter 2, two possible formulations for computing the force on a moving body, equa- 
tions (2.31) and (2.32). Wu [54] performed force calculations for a surface-piercing vessel 
using both equations and concluded that the two equations produce different answers. The 
correct expression for computations involving bounded fluids is equation (2.32). Equation 
(2.31) is appropriate if the body is surround by a fixed control surface on which U • n = 0 

[33]. 
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The time domain simulation provides the force time history in the same fashion as in 
the diffraction problem presented earlier in this thesis. The n th harmonic coefficients of the 
force are computed and related to the force at each order m using the method presented in 
Chapter 4. 
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Numerical Results 



The results for the radiation problem are compared with: 

• The body-exact linearized free surface formulation of Wu [52] for large amplitude 
oscillations; 

• The nonlinear BIEM of Silva and Peregrine [43]; 

• The linear analytic solution of Ogilvie (1963) for first order forces and the steady 
component of the second order force. 

Trajectory 

The cylinder is started from rest and the amplitude of the heave oscillation, £(£), is 
increased from zero to its limit state value using the following function [26]: 

= ( X - eX P _at ) Cax sin M) ( 5 - 4 ) 

where u is the frequency of oscillation and a is a startup constant. Equation (5.4) is used 
to provide a smooth startup. 

Radiation Condition- Numerical Beach 

The time simulations of the radiation problem are conducted until the force time history 
of the body reaches limit state. The method chosen to minimize the effect of transient 
behavior associated with the startup of the body from rest increases the time required for 
the problem to reach its limit state. Either a large computational domain or a method to 
eliminate wave reflection from the boundaries of the computational domain is required to 
ensure that the solution reaches its limit state. 
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The method developed by Dommermuth and Yue [8] and used by Liu [22] for enabling 
long time simulations for transient ship motions was used in this thesis. A tapering function 
was applied to the free surface elevation and free surface potential to smoothly truncate 
these quantities at the ends of the computational domain. The tapering function is shown 
in Figure 5-2. 




The tapering function is defined as [22]: 



1, A < x < L - A 

^(z.A) = { n((A-aO/A), 0 < £ < A 

U{(x-L + A) /A), L-A<x<L 

Where II (s) is the Hermitian polynomial [22] : 



(5.5) 
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n(s) = 1 - 462s 6 + 1980s 7 - 34655 s + 3080s 9 - 1386s 10 + 252s 11 , 0<s<l (5.6) 



The width of the tapering zone, A, is determined by physical and computational con- 
siderations. The tapering zone must be wider than the wavelength of the longest wave 
generated by the oscillating cylinder. If the tapering zone is shorter than this wavelength, 
the wave will be too strongly “reflected” from the domain boundaries. The difference be- 
tween the length of the computational domain and twice the length of the tapering zone, 
L* = L — 2 A, must be long enough to allow the simulation to reach and sustain physically 
relevant limit state behavior. Specifically, the ratio of the wave length of the wave gener- 
ated by the moving cylinder, A to L * must be small enough to allow multiple waves to exist 
within the domain. The wave length of the wave generated by the forced oscillatory motion 
at frequency u is estimated by the linear dispersion relation for gravity waves in a fluid of 
infinite depth, A = 2i Tg/u) 2 . 



Results 

The solutions for the free surface profile, the force time history, and the amplitudes of 
the steady and oscillatory forces generated by the combined HOSM-BIEM for the heaving 
circular cylinder are compared with the results generated from the other solution methods 
discussed in Chapter 1. 
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Figure 5-3 and Figure 5-4 show the free surface profile and force time history generated 
by Silva and Peregrine for a heaving circular cylinder. The parameters for the problem were 
(refer to Figure 5-1) : 



• R = 0.25; 

• H = -0.5; 

• d = 1.0; 

• f = 0.1; and, 



• Frequency of oscillation, u = 3.13 rad/sec. 




Figure 5-3: Free surface profile above heaving cylinder after approximately 5 periods of 
oscillation. The ordinate is the free surface elevation magnified by a factor of 500. The 
mean position of the cylinder center is (0.0, —0.5). (Silva and Peregrine, “Engineering 
Analysis with Boundary Elements”, 1990, Vol. 7, No. 4). 
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Figure 5-4: Time history of the vertical component of the force for the heaving cylinder. 
The ordinate is the normalized vertical force. The abscissa is time. (Silva and Peregrine, 
“Engineering Analysis with Boundary Elements”, 1990, Vol. 7, No. 4). 
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Figure 5-5 and Figure 5-6 show the free surface profile and force time history gener- 
ated by the combined HOSM-BIEM method for a heaving circular cylinder. The problem 
parameters were those listed on Page 83. 




Figure 5-5: Free surface profile above heaving cylinder after approximately 5 periods of 
oscillation. The ordinate is the free surface elevation magnified by a factor of 500. The 
mean position of the cylinder center is (L/ 2, —0.5). (Combined HOSM-BIEM method). 
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Figure 5-6: Time history of the vertical component of the force for the heaving cylinder. 
The ordinate is the normalized vertical force amplitude. The abscissa is time. (Combined 
HOSM-BIEM method). 
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The vertical scale for the free surface profiles presented in Figure 5-3 and Figure 5-5 
are magnified by a factor of 500 to facilitate the comparison of the results generated by the 
two different methods. The results generated by the BIEM method of Pergrine and Silva 
and the results generated by the present method are qualitatively similar. A harmonic 
analysis of the results of Peregrine and Silva was not performed. Therefore, a more detailed 
comparison of the results from the two methods could not be made. 

A comparison of the steady and first and second harmonics of the forces generated by 
the solution method of Wu [52] and the present method are presented in Figures 5-7, 5-8, 
and 5-9. Wu satisfied the body boundary condition on the exact position of the body and 
used linearized free surface dynamic and kinematic boundary conditions. Wu expressed the 
velocity potential in terms of a multipole expansion [52]: 



*s= £ £ 

m=l s=—oo 



o im0+isLJt 






+ £ £ K“" 



m— 1 $——oo 



o im0+isu)t 






(5.7) 



e im0+isut i r m j g the eX p ress i on f or the multipole expansion for a cylinder in an infinite fluid. 
/4i ( r , 9, t) and (r, 9, t) are required to enable the potential to satisfy the free surface 

boundary conditions [52]. 



87 




Figure 5-7: Vertical drift force on the heaving cylinder for kR = 0.1 ( ) and kR = 1.0 

( — ). Wu (1993) (□) ; present method (Q)- (H = 3 R). 




Figure 5-8: First harmonic of the vertical force on the heaving cylinder for kR = 0.1 ( ) 

and kR = 1.0 ( — ). Wu (1993) (□) ; present method (Q). ( H = 3/?). 
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Figure 5-9: Second harmonic of the vertical force on the heaving cylinder for kR = 0.1 
( ) and kR = 1.0 ( — ). Wu (1993) (□) ; present method (Q)- ( H = 3 R). 
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The results generated by the method of Wu and the present method agree well for 
heave amplitude to body radius ratios {£/R) less than 1.0. In the limit of vanishing heave 
amplitude to body radius ratios both methods converge to the linear solution of Ogilve [35]. 

The assumptions used by Wu to justify using the linear free surface boundary condition 
become invalid as the body approaches the free surface. For a deeply submerged body, the 
influence of the body on the free surface is of order 0(R/H) [52]. The linear free surface 
boundary condition ignores terms that involve products of order 0(R/H) 2 [52]. For heave 
amplitudes greater than 1.0 the body is close enough to the free surface to invalidate the 
deep submergence requirement for Wu’s method to be valid. In fact, for the real (physical) 
situation the waves break for the large orbital radii cases. 

5.2 Radiation Problem for a Circular Cylinder in a Circular 
Orbit 

The forces and free surface profiles resulting from the interaction of a fully submerged 
cylinder forced to undergo circular orbital motion with an initially quiescent free surface 
are presented in this section. The forces are evaluated as functions of: 

• The nondimensional orbital radius, 7 /i?; 

• The submergence ratio, h = H/R\ and, 

• The body radius to wavelength ratio (x 27r), kR = 2nR/X. 

The forces analyzed are: 

• The first and second harmonics of the vertical force, and 

• The steady vertical force. 
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The method for computing the forces on the body and the method for starting the 
problem from rest are those presented in the previous sections. 

Results 

The free surface profile, the force time history, and the amplitudes of the steady and 
oscillatory forces associated with the orbiting circular cylinder are compared with the results 
generated from the other solution methods discussed in Chapter 1. 

Figure 5-10 and Figure 5-11 show the free surface profile and force time history generated 
by Peregrine and Silva for the orbiting circular cylinder. The parameters for the problem 
were (refer to Figure 5-1) : 

• R = 0.3; 

• H = -0.5; 

• d = 1.0; 

• Orbital radius = 0.025; and, 

• Frequency of oscillation, u> = 3.13 rad/sec. 
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Figure 5-10: Free surface profile above the orbiting cylinder after approximately 5 periods 
of rotation. The ordinate is the free surface elevation magnified by a factor of 500. The 
cylinder center is orbiting about the point (0.0, —0.5). (Silva and Peregrine, “Engineering 
Analysis with Boundary Elements”, 1990, Vol. 7, No. 4). 
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Figure 5-11: Time history of the vertical ( ) and horizontal ( — ) components of the 

force for the orbiting cylinder. The ordinate is the normalized force amplitude. The abscissa 
is time. (Silva and Peregrine, “Engineering Analysis with Boundary Elements”, 1990, Vol. 
7, No. 4). 
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Figure 5-12 and Figure 5-13 show the free surface profile and force time history gen- 
erated by the new combined HOSM-BIEM method for the orbiting circular cylinder. The 
parameters for the problem were those listed on Page 91: 




Figure 5-12: Free surface profile above the orbiting cylinder after approximately 4 periods 
of oscillation. The ordinate is the free surface elevation magnified by a factor of 500. The 
cylinder center is orbiting about the point (L/2 , —0.5). (Combined HOSM-BIEM method). 
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Figure 5-13: Time history of the vertical ( ) and horizontal ( — ) components of the 

force for the orbiting cylinder. The ordinate is the normalized force amplitude. The abscissa 
is time. (Combined HOSM-BIEM method). 
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As was the case for the heaving cylinder, the comparison with the results of Peregrine 
and Silva provides a qualitative comparison of the results generated by these two different 
methods. 

It has been shown that at first order (c.f. Ogilvie [35]) and at second order (Wu [53]) 
an orbiting circular cylinder, in a fluid of infinite extent, generates waves in one direction; 
e.g., a cylinder undergoing a clockwise orbit generates only right going waves. Wu [52] 
demonstrated that if terms up to fifth order are included in a multipole expansion rep- 
resentation of the potential, the possiblity of waves being transmitted in both directions 
exist. However, these waves are small. The results of Silva and Peregrine and the present 
results demonstrate that in shallow water waves of significant amplitude are generated and 
transmitted in the “upstream” direction. 

This result is in contrast to that for a cylinder undergoing forced orbital motion in a 
fluid of infinite depth. Figure 5-14 shows the free surface profile above a circular cylinder 
orbiting beneath the free surface in a fluid of infinite depth. The problem parameters, with 
the exception of the fluid depth were those listed on Page 91. The cylinder generates waves 
in only the “downstream” direction when the fluid depth is infinite. 
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Figure 5-14: Free surface profile above orbiting cylinder in a fluid of infinite depth after 
approximately 5 periods of oscillation. The ordinate is the free surface elevation magnified 
by a factor of 500. The cylinder is orbiting about the point (L / 2 , —0.5). (Combined 
HOSM-BIEM method). 
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A comparison of the steady and first and second harmonics of the forces generated by 
the solution method of Wu [52] and the present method are presented in Figures 5-15, 5-16, 
and 5-17. 




Figure 5-15: Vertical ( ) and horizontal ( — ) drift forces on orbiting circular cylinder. 

Wu (1993) (□); present method (O). (H = 3 R, kR = 0.5). 

As was the case for the heaving cylinder the two methods agree well for small amplitude 
of oscillation (orbital radius) to cylinder radius ratios. 

For large orbital radii, the cylinder generates waves with local steepnesses that exceed 
the limiting Stokes steepness. Figure 5-18 shows the local wave steepness for a case where 
the orbital radius to body radius is 1.5 and the mean depth is 3 times the body radius. For 
large orbital radii, the method becomes unstable as the body motion generates excessively 
steep waves, and does not reach limit state behavior before the simulation stops. 
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Figure 5-16: First harmonic of the vertical ( ) and horizontal ( — ) forces on orbiting 

circular cylinder. Wu (1993) (□); present method (Q)- (H = 3 R, kR = 0.5). 




Figure 5-17: Second harmonic of the vertical ( ) and horizontal ( — ) forces on orbiting 

circular cylinder. Wu (1993) (□); present method (O)- (H = 3 R, kR = 0.5). 
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Figure 5-18: Free surface profile above the orbiting cylinder. (H = 3 R\ kR = 0.5). 
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5.3 Summary 



The ability of the combined high-order spectral- boundary integral equation method to sim- 
ulate large amplitude oscillatory motions of a circular cylinder beneath a free surface has 
been demonstrated. The method enables forces and free surface profiles to be evaluated up 
to limiting Stokes’ wave steepness for a range of problems of interest. Problems involving 
either finite or infinite depth can be modelled. The importance of using a solution method 
which identifies the existence of physically impossible waves is demonstrated by the applica- 
tion of the present method to large amplitude body motions near a free surface. A method 
such as the one used by Wu, which does not account for non-linear free surface interactions 
and locally steep waves, extends linear results into flow regimes where linear assumptions 
are invalid. 

In the next chapter, the HOSM-BIEM method will be applied to the combined radiation- 
diffraction problem for a fully submerged cylinder free to respond to incident waves. 
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Chapter 6 



Combined Radiation and 
Diffraction Problem for a 
Submerged Circular Cylinder 



The purpose of this chapter is to demonstrate the application of the combined high order 
spectral-boundary integral equation method to the problem of a fully submerged circular 
cylinder free to respond to an incident wave field. Figure 6-1 defines the problem parameters. 



102 




Figure 6-1: Problem parameters for the combined problem 



• R = the body’s radius 

• r 0 = the radius of body’s orbital trajectory 

• (x 0 , h 0 ) = the location of the point about which the body is orbitting 

• A = the wavelength of the ambient waves 

• A = the wave amplitude 



103 



6.1 The Combined Radiation and Diffraction Problem for a 



Circular Cylinder Free to Respond to Waves 

The analyses of neutrally bouyant circular cylinders free to respond to incident wave fields 
were facilitated by the generation of cylinder trajectories. The cylinder trajectories were 
generated by: 

• Specifying initial conditions for the cylinder position, velocity and acceleration; 

• Solving for the instantaneous force on the cylinder using the solution methodology 
detailed in Chapters 2 and 3; and, 

• Updating the cylinder position, velocity, and acceleration using the state varible ap- 
proach represented by Equations 2.33, 2.34, 2.35, 2.36, and 2.37. 

Initial Conditions 

Ogilvie [35] used linear theory to show that the amplitude and phase angle of the orbital 
motion of a submerged circular cylinder differs from that of a water particle centered at 
the location of the cylinder’s origin by a quantity of fourth order in kR, where k is wave 
number and R is the body’s radius. This fact was used to estimate an intitial velocity and 
accleration for the cylinder. 

The first-order orbital radius of a water particle located initially at ( x 0) h 0 ) (see Figure 6- 
1) would be: 



r o (0) = Ae~ kh °e~ iwt \t=o = Ae~ kho 



where to is the wave frequency. 



( 6 . 1 ) 
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The initial velocity {v 0 (t = 0) = (u(0),u>(0))) and acceleration (a 0 (t = 0) = (a x (0), 0^(0))) 
of the cylinder were set the same as that of a fluid particle at the cylinder center. 



u(0) = ujAe~ ktl ° cos (kx) 


(6.2) 


w(0) = uAe~ kh ° sin(A:a:) 


(6.3) 


a x (0) = u 2 Ae~ kh ° sin(kx) 


(6.4) 


a y ( 0) = —uj 2 Ae~ ktl ° cos (kx) 


(6.5) 



Method of Computing Forces and Generating the Trajectory of the Cylinder 

The potential on the body is computed at each time step by solving the initial-boundary 
value problem presented in Chapters 2 and 3. The potential is used to compute the instan- 
taneous force on the body, F(t), by the application of Equation 5.3. The computed force 
is used to update the body’s acceleration, velocity, and position using the state variable 
representation of these quantities presented in Chapter 2: 



a„(t) 


= X(t) = F(t)/M 


(6.6) 








v 0 (t + St) 


= V 0 (t) + J a 0 {t)dt 


(6.7) 




rt+5t 




a 0 [t + St) 


= a 0 (t) + j v 0 {t)dt 


(6.8) 



Numerical Results 

The trajectories of submerged cylinders beneath waves were simulated. The waves used 
in the simulations were waves travelling in the positive X direction (left to right in the 
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figures that follow). 

The convergence with respect to time step size on the trajectory of a cylinder free to 
respond to waves is shown in Figure 6-2. Figure 6-2 shows that in order to achieve accurate 
cylinder trajectories a time step size of T/ 256 or smaller is required. 




Figure 6-2: Convergence of the trajectory of a submerged circular cylinder with respect to 
time step size, dt. Wave slope kA = 0.08; cylinder initial position ( x o ,h 0 ) = (157T, — 1.2); 
kR = 0.4. The trajectory of the cylinder center is shown. Order M = 4, number of Fourier 
modes N x = 2048, domain size Nw = 16, Nfu ter = 24 Nw, Nb = 128. 
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The convergence with respect to number of body segments on the trajectory of a cylinder 



free to respond to waves is shown in Figure 6-3. 




47.08 47.10 47.12 47.14 47.16 47.18 

X 



Figure 6-3: Convergence of the trajectory of a submerged circular cylinder with respect to 
number of body segments, Nq . Wave slope kA = 0.08; cylinder initial position (x 0 ,/i 0 ) = 
(157r, — 1.2); kR = 0.4. The trajectory of the cylinder center is shown. Order M = 4, 
number of Fourier modes N x = 2048, domain size Nw = 16, Nju ter = 24 Nw, dt = T/ 256. 
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Figures 6-4 and 6-6 show the trajectories for two different cylinders free to respond to 
incident waves. The body radius for the problem represented by Figure 6-4 is 2.5 times the 
body radius for the problem represented by Figure 6-6. The expanded cylnder trajectories 
for each case are shown in Figures 6-5 and 6-7. Figures 6-4 and 6-6 show the cylinder 
positions at three different time steps. The free surface profiles shown are the free surface 
profiles above the cylinder when the cylinders’ positions are as shown by the solid lines. 




Figure 6-4: Cylinder trajectory and an associated free surface profile. Cylinder position at 

t = 0 ( ); cylinder position at t = 3T ( — ); cylinder position at t = 3.75 T (— • — • — ). 

Wave slope kA = 0.08; kR = 1.0; cylinder initial position (x 0 , h 0 ) = (L/2, — 1.25). Order 
M = 4, number of Fourier modes N x = 2048, domain size Nw = 16, Nfut er = 24 Nw, 
dt = T/ 256, N b = 192. 
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Figure 6-5: Expanded cylinder trajectory. Wave slope kA = 0.08; kR = 1.0; cylinder initial 
position {x 0 ,h Q ) = (L/ 2,— 1.25). Order M = 4, number of Fourier modes N x = 2048, 
domain size Nw = 16, Nfut er = 24 Nw, dt = T/ 256, Ng = 192. 
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Figure 6-6: Cylinder trajectory and an associated free surface profile. Cylinder position at 

t = 0 ( ); cylinder position at t = 3T (— • — • — ); cylinder position at t = IT. Wave 

slope kA = 0.08; kR = 0.4; cylinder initial position (x 0 ,h 0 ) = (L/2, —1.25). Order M = 4, 
number of Fourier modes N x = 2048, domain size Nw = 16, Nfu ter = 24 Nw, dt = Tj 256, 
N b = 192. 
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Figure 6-7: Expanded cylinder trajectory. Wave slope kA = 0.08; kR = 0.4; cylinder initial 
position (x 0 , h Q ) — (L/ 2, —0.8). Order M = 4, number of Fourier modes N x = 2048, domain 
size ATw = 16, Nfm er = 24A/V, dt = Tj 256, AT# = 128. 
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The cylinder trajectories reflect the influence of the imposed initial conditions in that 
for both cases the cylinders initially follow the trajectory that would be assumed by a fluid 
particle located at the cylinder’s initial position (x 0 ,y 0 ). As time increases the steady and 
oscillatory forces induced by the fluid-body interaction generate motions which differ from 
that of a fluid particle. For the large body the presence of a negative horizontal drift motion 
is apparent while for the smaller body this does not appear to be the case. 

The simulation of the interaction between the large cylinder and the free surface was 
terminated due to the development of short, locally steep waves as the body became close 
to the free surface. The presence of the small locally steep waves is apparent in Figure 6-4. 
The development of locally steep waves helps to define the limitations of the method of 
this thesis and to identify flow regimes where linear theories used for modelling wave-body 
interactions are invalid. The free surface boundary conditions used in this thesis can not 
be used to simulate breaking waves. 

6.2 Summary 

The combined high-order spectral-boundary integral equation method has been successfully 
applied to the problem of a fully submerged neutrally bouyant cylinder free to respond to 
incident waves. The method is stable enough to enable simulations of sufficient duration to 
evaluate the steady and oscillatory responses of interest. 
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Chapter 7 



Conclusions and Recommendations 



A combined high-order spectral method boundary integral equation method was developed 
for the analysis of the nonlinear interaction of non-breaking waves with fully submerged 
bodies of arbitrary geometry undergoing arbitrary motions. The motivation for the devel- 
opment of the method wa s the desire to improve the computation speed for wave interaction 
problems involving bodies of arbitrary geometry undergoing arbitrary motions near the free 
surface. 

The method developed in this thesis represents the total fluid potential as an expansion 
in a perturbation series of terms proportional to wave slope. The total potential is the sum 
of a spectral potential and a body potential. The method uses a spectral representation of 
free surface quantities and a boundary element representation of the body. The spectral 
representation of the free surface elevation and potential enables the fast computation of 
free surface quantities. The boundary element representation of the body enables bodies of 
arbitrary shape to be modelled. Fully nonlinear free surface boundary conditions are used 
as time evolution equations for the free surface elevation and free surface potential. 

The effectiveness of the method was demonstrated through the study of three classes 
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of problems: diffraction of incident waves by a stationary body, radiation from a body 
undergoing forced oscillatory motion beneath a free surface, and a combined radiation- 
diffraction problem for a body free to respond to an incident wave field. 

The success and limitations of the present method in computing the forces associated 
with nonlinear wave body interactions are demonstrated in Chapters 4, 5 , and 6. The 
present method succeeds in computing 

• Oscillatory forces; 

• Steady drift forces; 

• Free surface profiles up to limiting wave steepness; and, 

• Trajectories for bodies free to respond to waves. 

The forces and free surface profiles produced by the present method for the diffraction 
and radiation problems compare well with those produced by other methods restricted to 
interactions involving non-breaking waves, and, for waves of small steepness, compare well 
with those produced by linear methods. Through evaluation of the free surface profiles 
generated by the wave-body interaction, locally steep profiles which exceed limiting wave 
steepness can be identified and used to determine when methods based on perturbation 
expansions are no longer appropriate or accurate for the characterization the posed physical 
problems. 

The ability of the present method to evaluate bodies of arbitrary geometry was demon- 
strated by the evaluation of a body of pontoon cross section. The ability of the present 
method to evaluate bodies undergoing arbitrary oscillatory motions was demonstrated by 
evaluating the large amplitude oscillatory motions of a circular cylinder. 
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Recommendations 



The present method is formulated to accomodate non-breaking waves. As shown in the 
analysis of the diffraction and radiation problems, locally steep waves limit the ability of 
the present method to evaluate wave-body interactions when the body is close to the free 
surface even when the waves in the ambient wave field are not steep. 

The periodic nature of the formulation of the present method limits problem simulation 
time. For the radiation problem the use of a numerical damping zone enables longer simlu- 
ation times. A method for enabling long time simulations for the diffraction and combined 
problems needs to be developed to enable the method to perform long time simulations. 
This is important for the manuevering simulations where the distance travelled by a body 
is more than a few body lengths. 

The present method was demonstrated for two dimensional non-lifting bodies. The 
extension of the method to three dimensions and problems involving lift is in principle 
straightforward, but remains to accomplished and demonstrated. 
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